#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<gmp.h>
#include<sys/time.h>
#include<pthread.h>

// Definições
#define NUMDIGITOS 10000001
#define NOMEARQUIVO "gauss_legendre_paralelo.txt"

// Definições auxiliares
#define LOG_2_10 3.32192809488736234787

// Estrutura para o timer
struct timeval timerInicial, timerFinal;

// Variáveis do método
mpf_t a, b, t, tSub, p;
mpf_t aN, bN, tN, pN;
mpf_t pi;

// Calcula A(n+1)
void *calculaAn(void *args) {
    // Calcula
    mpf_add(aN, a, b);
    mpf_div_ui(aN, aN, 2);

    // Retorno
    return NULL;
}

// Calcula B(n+1)
void *calculaBn(void *args) {
    // Calcula
    mpf_mul(bN, a, b);
    mpf_sqrt(bN, bN);

    // Retorno
    return NULL;
}

// Calcula T(n+1)
void *calculaTn(void *args) {
    // Calcula
    mpf_sub(tSub, a, aN);
    mpf_pow_ui(tSub, tSub, 2);
    mpf_mul(tSub, tSub, p);
    mpf_sub(tN, t, tSub);

    // Retorno
    return NULL;
}

// Calcula P(n+1)
void *calculaPn(void *args) {
    // Calcula
    mpf_mul_ui(pN, p, 2);

    // Retorno
    return NULL;
}

// Função principal
int main(void) {
    // Timer inicial
    gettimeofday(&timerInicial, NULL);

    // Variáveis de iteração
    int i, numIteracoes;

    // Threads
    pthread_t threadAn, threadBn, threadTn, threadPn;

    // Número de iterações necessárias = log2(númeroDeCasasDesejadas)
    numIteracoes = log(NUMDIGITOS) / log(2.0);

    // Precisão
    mpf_set_default_prec(LOG_2_10 * NUMDIGITOS);

    // Inicializa as variáveis
    mpf_inits(a, b, t, p, pi, NULL);
    mpf_inits(aN, bN, tN, pN, tSub, NULL);

    // Arquivo de saída
    FILE *arqSaida;
    arqSaida = fopen(NOMEARQUIVO, "w");

    // ERRO: Arquivo inexistente
    if(arqSaida == NULL)
        exit(EXIT_FAILURE);

    // MÉTODO DE GAUSS-LEGENDRE
    // Condições iniciais
    // A(0) = 1
    mpf_set_d(a, 1.0);

    // B(0) = 1 / sqrt(2)
    mpf_set_d(b, 2.0);
    mpf_sqrt(b, b);
    mpf_ui_div(b, 1, b);

    // T(0) = 1/4
    mpf_set_d(t, 0.25);

    // P(0) = 1
    mpf_set_d(p, 1.0);

    // Imprime o cabeçalho
    printf("VALOR DE PI\n\n");

    // Método
    for(i = 0; i < numIteracoes; i++) {
        // A(n+1) = A(n) + B(n) / 2
        pthread_create(&threadAn, NULL, calculaAn, NULL);

        // B(n+1) = sqrt( A(n) * B(n) )
        pthread_create(&threadBn, NULL, calculaBn, NULL);

        // P(n+1) = 2 * P(n)
        pthread_create(&threadPn, NULL, calculaPn, NULL);

        // T(n+1) = T(n) - P(n) * ( A(n) - A(n+1) ) ^ 2
        // T(n+1) depende de A(n+1)
        pthread_join(threadAn, NULL);
        pthread_create(&threadTn, NULL, calculaTn, NULL);
        pthread_join(threadTn, NULL);

        // PI = ( (A(n) + B(n)) ^ 2 ) / (4 * t(n))
        mpf_add(pi, a, b);
        mpf_pow_ui(pi, pi, 2);
        mpf_div_ui(pi, pi, 4);
        mpf_div(pi, pi, t);

        // Imprime o valor aproximado de PI
        gmp_printf("PI(%d) = %.15Ff\n", i, pi);

        // A(n) = A(n+1)
        // A(n) depende de A(n+1)
        mpf_set(a, aN);

        // B(n) = B(n+1)
        // B(n) depende de B(n+1)
        pthread_join(threadBn, NULL);
        mpf_set(b, bN);

        // T(n) = T(n+1)
        // T(n) depende de T(n+1)
        mpf_set(t, tN);

        // P(n) = P(n+1)
        // P(n) depende de P(n+1)
        pthread_join(threadPn, NULL);
        mpf_set(p, pN);
    }

    // Salva o pi
    mpf_out_str(arqSaida, 10, NUMDIGITOS, pi);

    // Limpa as variáveis
    mpf_clears(a, b, t, p, pi, NULL);
    mpf_clears(tSub, aN, bN, tN, NULL);

    // Encerra o arquivo
    fclose(arqSaida);

    // Timer final
    gettimeofday(&timerFinal, NULL);

    // Imprime o resumo
    printf("\n\nRESUMO DA EXECUCAO\n\n");

    // Imprime o número de iterações
    printf("Numero de iteracoes: %d\n", numIteracoes);

    // Imprime o tempo de execução em microssegundos (us)
    printf("Tempo de execucao: %d us\n", (int) (timerFinal.tv_usec - timerInicial.tv_usec));

    // Retorno
    return EXIT_SUCCESS;
}
